//  Copyright (c) 2012 M.A. (Thijs) van den Berg, http://sitmo.com/
//
//  Use, modification and distribution are subject to the MIT Software License. 
// (See accompanying file LICENSE.txt or copy at http://www.stdfin.org/LICENSE.txt)
// 

#ifndef STDFIN_PRICING_GBM_OPTION_VANILLA_HPP
#define STDFIN_PRICING_GBM_OPTION_VANILLA_HPP

#include <cmath>
#include "pricing_math.hpp"

namespace stdfin{

/**
 * Calculate the price of a vanilla call option.
 *
 * A vanilla call option has a pay-off max(S-X,0) at expiration.
 *
 * @param S
 *   Current value of the underlying.
 * @param Y
 *   Yield of the underlying.
 * @param v
 *   Annualized volatility of the underlying.
 * @param r
 *   Continuous compounded interest rate.
 * @param X
 *   Strike of the vanilla call option
 * @param t
 *   Time in years till expiration.
 *
 * @return
 *   The price of the call option.
 */
inline double gbm_vanilla_option_call_price(double S,double Y,double v,double r,double X,double t)
{
    double d1 = (  log(S/X) + (Y + 0.5*v*v)*t  ) / (v*sqrt(t));
    double d2 = d1 - v*sqrt(t);
    return S*exp((Y-r)*t)*pricing::N(d1) - X*exp(-r*t)*pricing::N(d2);
}

/**
 * Calculate the price of a vanilla put option.
 *
 * A vanilla put option has a pay-off max(X-S,0) at expiration.
 *
 * @param S
 *   Current value of the underlying.
 * @param Y
 *   Yield of the underlying.
 * @param v
 *   Annualized volatility of the underlying.
 * @param r
 *   Continuous compounded interest rate.
 * @param X
 *   Strike of the vanilla call option
 * @param t
 *   Time in years till expiration.
 *
 * @return
 *   The price of the put option.
 */
 inline double gbm_vanilla_option_put_price(double S,double Y,double v,double r,double X,double t)
{
    double d1 = (  log(S/X) + (Y + 0.5*v*v)*t  ) / (v*sqrt(t));
    double d2 = d1 - v*sqrt(t);
    return -S*exp((Y-r)*t)*pricing::N(-d1) + X*exp(-r*t)*pricing::N(-d2);
}

inline double gbm_vanilla_option_call_delta(double S,double Y,double v,double r,double X,double t)
{
    double d1 = (  log(S/X) + (Y + 0.5*v*v)*t  ) / (v*sqrt(t));
    return exp((Y-r)*t)*pricing::N(d1);
}

inline double gbm_vanilla_option_put_delta(double S,double Y,double v,double r,double X,double t)
{
    double d1 = (  log(S/X) + (Y + 0.5*v*v)*t  ) / (v*sqrt(t));
    return -exp((Y-r)*t)*pricing::N(-d1);
}

inline double gbm_vanilla_option_gamma(double S,double Y,double v,double r,double X,double t)
{
    double d1 = (  log(S/X) + (Y + 0.5*v*v)*t  ) / (v*sqrt(t));
    return exp((Y-r)*t)*pricing::n(d1) /   ( S*v*sqrt(t) );
}

inline double gbm_vanilla_option_call_theta(double S,double Y,double v,double r,double X,double t)
{
    double d1 = (  log(S/X) + (Y + 0.5*v*v)*t  ) / (v*sqrt(t));
    double d2 = d1 - v*sqrt(t);
    return  (Y-r) * exp((Y-r)*t) * S * pricing::N(d1)
               +r * exp(-r   *t) * X * pricing::N(d2)
                  + exp((Y-r)*t) * S * v * pricing::n(d1) / (2*sqrt(t));
}

inline double gbm_vanilla_option_put_theta(double S,double Y,double v,double r,double X,double t)
{
    double d1 = (  log(S/X) + (Y + 0.5*v*v)*t  ) / (v*sqrt(t));
    double d2 = d1 - v*sqrt(t);
    return  (r-Y) * exp((Y-r)*t) * S * pricing::N(-d1)
           -r * exp(-r   *t) * X * pricing::N(-d2)
              + exp((Y-r)*t) * S * v * pricing::n(d1) / (2*sqrt(t));
}

inline double gbm_vanilla_option_vega(double S,double Y,double v,double r,double X,double t)
{
    double d1 = (  log(S/X) + (Y + 0.5*v*v)*t  ) / (v*sqrt(t));
    double d2 = d1 - v*sqrt(t);
    return exp((Y-r)*t) * S * sqrt(t) * pricing::n(d1);
}

inline double gbm_vanilla_option_call_rho(double S,double Y,double v,double r,double X,double t)
{
    double d1 = (  log(S/X) + (Y + 0.5*v*v)*t  ) / (v*sqrt(t));
    double d2 = d1 - v*sqrt(t);
    return t * exp(-r*t) * X * pricing::N(d2) - t * exp((Y-r)*t) * S * pricing::N(d1);
}

inline double gbm_vanilla_option_put_rho(double S,double Y,double v,double r,double X,double t)
{
    double d1 = (  log(S/X) + (Y + 0.5*v*v)*t  ) / (v*sqrt(t));
    double d2 = d1 - v*sqrt(t);
    return -t * exp(-r*t) * X * pricing::N(-d2) + t * exp((Y-r)*t) * S * pricing::N(-d1);
}

inline double gbm_vanilla_option_call_yield(double S,double Y,double v,double r,double X,double t)
{
    double d1 = (  log(S/X) + (Y + 0.5*v*v)*t  ) / (v*sqrt(t));
    double d2 = d1 - v*sqrt(t);
    return t * exp((Y-r)*t) * S * pricing::N(d1);
}

inline double gbm_vanilla_option_put_yield(double S,double Y,double v,double r,double X,double t)
{
    double d1 = (  log(S/X) + (Y + 0.5*v*v)*t  ) / (v*sqrt(t));
    double d2 = d1 - v*sqrt(t);
    return -t * exp((Y-r)*t) * S * pricing::N(-d1);
}

} // namespace

#endif